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Abstract 

This paper introduces a class of approximate transparent boundary conditions 
for the solution of Helmholtz-type resonance and scattering problems on un- 
bounded domains. The computational domain is assumed to be a polygon. A 
detailed description of two variants of the Hardy space infinite element method 
which relays on the pole condition is given. The method can treat waveguide- 
type inhomogeneities in the domain with non-compact support. The results of 
the Hardy space infinite element method are compared to a perfectly matched 
layer method. Numerical experiments indicate that the approximation error of 
the Hardy space decays exponentially in the number of Hardy space modes. 
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1. Introduction 

To solve numerically the Helmholtz equation 

- Au(x) - k 2 n(x)u(x) = 0, x E tt (1) 

on an unbounded domain fi with some boundary conditions on d£l and a radia- 
tion condition at infinity, the computational domain is typically restricted to a 
bounded interior f2i n t. Here n(x) is the refraction index and k is the wavenum- 
ber. In case of a scattering problem k > is given, whereas in case of a resonance 
problem k with positive real part is the sought resonance. Applying transparent 
boundary conditions at the artificial interface T of the exterior/interior domain, 
the problem can be restricted to solving the Hclmholtz equation on the interior 
domain only. Transparent boundary conditions on T have to model the correct 
radiation condition at infinity. 
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For homogeneous exterior domains and k > the correct radiation condition 
is the Sommcrfeld radiation condition. It is known that the exact transparent 
boundary condition, i.e. the Calderon or Dirichlet-to-Neumann map, is non- 
local. Currently used methods to approximate or realize transparent boundary 
conditions are based on separable coordinates and special functions [7J, infinite 
elements [TJH], perfectly matched layer (PML) constructions [3^112113], boundary 
integral approaches [14] and local high order approximations [5]. 

Except for some PML methods they depend non- linearly on k 2 , which is a 
severe drawback when solving resonance problems, where one seeks non-trivial 
eigenpairs (it, n 2 ) for ([!]) with vanishing boundary conditions on 9f2 and a ra- 
diation condition at infinity. Using finite elements to discretize ([I]) in i? 1 (J7; n t) 
leads to a generalized eigenvalue problem. Employing transparent boundary 
conditions, that are non-linear in the eigenvalue, would make the eigenvalue 
problem non-linear. Although it is possible to solve the resulting non-linear 
eigenvalue problems, see e.g. [IS], it is reasonable to avoid them. Therefore 
PML methods are currently the standard method for solving resonance prob- 
lems, see e.g. [5J [T7J. Under the name complex scaling they have been used 
since the 1970s for the theoretical study and the numerical computation of res- 
onances in molecular physics |26j . Unfortunately these methods give rise to 
spurious resonance modes and several parameters have to be optimized for each 
problem. The PML method we use to check the results of our Hardy space 
infinite element method choses the thickness and the discretization of the layer 
adaptively [27] 132]. 

The theoretical framework of the Hardy space infinite element (HSIE) method 
is the pole condition by F. Schmidt [23] [24]. The pole condition considers the 
Laplacetransform of the exterior solution with respect to some generalized dis- 
tance variable. A solution is then called (purely) outgoing if this Laplacetrans- 
form has no singularity in the lower complex half plane, vice versa a solution is 
(purely) incoming if its Laplacetransform has no singularity in the upper com- 
plex half plane. In [13] it is shown that for homogeneous exterior domains this 
condition on the singularities of the Laplacetransform, i.e. all its singularities 
are located in the upper complex half, is equivalent to the Sommerfeld radiation 
condition. 

Compared to former numerical realizations of the real axis approach of the 
pole condition [12 , which were based on BDF and Runge-Kutta methods, the 
HSIE method which is based on a Galerkin method in the Hardy space H + (D) 
of the complex unit disk, shows exponential convergence and is simple to imple- 
ment in finite element codes. It was first presented in [2UJ [TT] for homogeneous 
exterior domains with spherical interface T. The cut function approach of the 
pole condition |12j was discretized in |25j using a collocation method in some 
sort of Hardy space and shows exponential convergence in experiments, too. 
Moreover it allows for the evaluation of the exterior field. However this ap- 
proach is not linear in k 2 and exhibits problems concerning the stability under 
perturbations of the boundary data. 

Extending [TT] we assume here that the interface T is the boundary of a 
convex polygon P, i.e. Q- m t = P D f2 and f2 cxt = M 2 \ P, and that the exte- 
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rior domain f2 oxt is discretized by infinite trapezoids, such that the refractive 
index n(x) of ([!]) is constant on each trapezoid. These trapezoids are the im- 
ages of bilinear mappings from a reference strip. The infinite direction of this 
reference strip is mapped into the Hardy space H + (D), where a L 2 -orthogonal 
basis is given by trigonometric monomials. The basis functions of our Galerkin 
method are therefore tensor products of standard finite element functions with 
the trigonometric monomials in the Hardy space. 

The method presented here treats scattering problems as well as the corre- 
sponding resonance problems. The pole condition as a mean to realize trans- 
parent boundary conditions for time dependent problems is considered in |21j . 
The discretization used there is almost equivalent to one employed here. 

In section 2 we present the HSIE method from a practical point of view first 
in one dimension and then in two dimensions. In section 3 we shortly review 
the PML method used here. Numerical results are presented in section 4 with 
a comparison of both methods. 



2. Hardy space infinite element method 

To explain the basic ideas of HSIEs we shortly discuss one-dimensional 
problems, even though for such problems there exist simpler and more effi- 
cient methods to treat the unboundedness of the domain. Nevertheless, as the 
multi-dimensional elements are tensor products of standard finite elements and 
one-dimensional infinite elements, it is useful to start with the simple case. 

2.1. One- dimensional elements 

We consider the Hclmholtz equation 

- u"(r) - K 2 n(r)u(r) = 0, r > 0, (2a) 
u'(0) = g, (2b) 
u is outgoing (2c) 

with complex wave number k€C with positive real part (K(/t) > 0), boundary 
value g e C, and positive potential n € L°°((0, oo)) satisfying n(r) = 1 for r > a 
for some a. The Sommerfeld radiation condition 



lim [d r u{r) — inu(r)) = 

r— >oo 

guarantees for k with nonnegative imaginary part (3(k) > that is well- 
posed and that the solution u is outward radiating. Solutions to (2a) may be 
decomposed into an interior part u; n t := u|[o o] an d & n exterior part u cxt (r) := 
u(r + a), r > 0, with 

Uext (r) = C x e iKT + C 2 e- lKr and d+C 2 = u int (a). 

The term C\e %KT corresponds to an outgoing radiating wave, that satisfies the 
Sommerfeld radiation condition, and Cie~ lKr corresponds to an incoming wave. 
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Therefore C2 = and hence u e xt(f) = Ui nt (a)e lKr . The solution of ^ restricted 
to [0, a] is thus given by the simple boundary value problem 

-Ui nt (r)-n(r)K 2 u int (r)=0, u[ nt (0) = g, u[ nt (a) = iKu int (a). (3) 

Thus in the one-dimensional case the Sommerfeld radiation condition is satisfied 
at the boundary of the interior domain. Therefore it can be used for k with 
negative imaginary part, too. 

Note that the resonance problem corresponding to ([3| leads to a quadratic 
eigenvalue problem in n whereas the discretization using Hardy space infinite 
elements will lead to an eigenvalue problem that is linear in k 2 . 

Hardy space infinite elements rest on the fact, that the Laplace transform 
of the exterior solution 

poo 

(Cu cxt )(s) := / e~ sr u C xt(r)dr, > |3(«)| 

Jo 

has a holomorphic extension except for two poles at ±in 

C \C ie lK ' + C 2 e- lK '\ (s) = + - C ^. 

s — in s + in 

Hence, u ox t is outgoing if and only if £it ex t has no poles with negative imaginary 
part. An equivalent formulation in terms of an appropriate Hardy space is given 
in Definition [5] 

Definition 1 (Hardy space). Let P~ o = {s e C : ^s(s/kq) < 0} be the half 
plane below the line KqR. through the origin and k , see Fig. |2.1| The Hardy 
space H~(P~ o ) is the space of all functions /, that are holomorphic in P K(1 , such 
that 

\f(K x - K ie)\ 2 dx 



is uniformly bounded for e > 0. 

Let D = {s G C : \s\ < 1} be the open unit disk. The Hardy space H + (D) is 
the space of all functions /, that are holomorphic in D, such that 

^ \f(re^)\ 2 dt 



is bounded uniformly for re [0, 1). 

Due to the uniform boundedness of / there exist in both cases a L 2 function 
on the boundary, which is uniquely determined by / and which determines vice 
versa uniquely the function in the domain. Hence, we identify a Hardy space 
function / £ H~(P~ o ) or / G H + {D) with its boundary function / e L 2 (kqM.) 
and / e L 2 (S 1 ) respectively. 

Definition 2 (Pole condition). Let kq be a complex constant with positive 



real part and 9?(k/ko) > 0. Then a solution u to (2a I is said to obey the pole 
condition and is called outgoing, if the holomorphic extension of the Laplace 
transform of the exterior part lies in the Hardy space H~{P~ ). 
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p- 



Figurc 1: sketch to definition [T] 



The constant kq will act as a tuning parameter in the method presented 
below. In future we omit the formulation "holomorphic extension of the Laplace 
transform" and shortly write Laplace transform C. The Hardy space H~{P~ Q ) 
is a Hilbert space [5J [TU] , with the standard L 2 -norm, and the following Lemma 
connects the two Hardy spaces H~(P~ ) and H + (D). 

Lemma 3 (Mobius transform). The mapping 

M K0 : H-(P- )^H+(D) : / ^ (M K J)(z) := / (i«o^l > ) (4) 
is up to a factor ^/2|kq| unitary. 



z-l z-l 



Due to the explicit knowledge of u ex t the transformed function U :— M Ko Cu cxt 
is given by 

m = . . + u ;f a) ( n = ± ( K —^Y *>. ( 5 ) 

IK0(Z + 1) - IK{Z - 1) 1{K + KQ) \K+ K J 



Since 



K— Kp 



< 1, we could expect exponential convergence for the exterior 

solution, if we use the first N + 1 trigonometric monomials {z° , z 1 , z N } as a 
Galerkin basis of the space H + {D). For the interior part iti n t £ i? 1 ([0, a]) we 
use a standard finite element method. 

Both methods have the term Ui n %(a) in common and we call the associated 
boundary degree of freedom uq. As usual it couples the finite elements for the 
interior part with the infinite elements for the exterior part. Note, that in the 
formulation (J5| all degrees of freedom for U would couple with Uj n t(a), since 
Uint(a) = 2iKolI(l). In order to get a local coupling of the boundary degree of 
freedom and the inner degrees of freedom, we decompose U = j^T-(uq,U) t 
with 

T-(^yz)--=\{uv + {z-l)U{z)) (6) 
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for (uo, U) T G C x H + (D) and use the trigonometric monomials as a basis for 
U. 

It remains to derive a variational formulation for (u- m t,U) T in ii^QOja]) x 
H + (D), which was done in [TT]. Here i? 1 ([0, a]) is the Sobolev space of weakly 
differentiable functions on the interval [0, a]. The basic idea is to use the prop- 
erties of the Fourier transform to get the identity 



f(r)g(r)dr = -2i Ko A(M Ko Cf,M Ko Cg) 



with 



A(F,G) 



2n 



S 1 



F(z)G(z)\dz\, 



F,G&H + (D). 



(7) 



(8) 



This holds for u cxt and suitable test functions v cxt , as well as for the deriva- 
tives u' cxt and v' cxt . Using the decomposition A4 Ko £u ex t — ik^7~-(uo,U) T and 
M. Ko Cv cxt = jj^7~-{vq, V) T we obtain simple formulas for the derivatives of u ex t 
and Uext 

M K0 Cf'=T+ ( f °) withT+ 



(*) == 2 (/o + (* + !)*■(*)) • (9) 



Now we are able to deduce from the formal variational formulation of ( 2a | 



(u' int v' int - K 2 nu int v in t)dr + / (u' cxt v' cxt - n 2 nu cxt v cxt )dr = -gv int (0) 



and ([7]) the variational equation in i? 1 ([0, a]) x H + (D): 



with 



B 



Mint 
£7 



2ikqA T-i 



B 



fint 

r 



Mint \ /" «int 

17 H V 



-flWint(O) 



(10) 



:= / (u int W int - K Mint w int )dr 



,71 



'Ml 



-A 71 



Since for the trigonometric monomials A^ , z k ^j — Sj^, the implementation of 
the exterior part of the bilinear form B reduces to the implementation of the 
operators 7± : C x H + (D) — > H + (D), if the finite dimensional ansatz space 
IIjv := span-jl , z 1 , z N } is used for H + (D): 



( 1 ±1 



N,± 



\ 

1 ±1 

1 / 



(11) 
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(a) (b) 

Figure 2: a) 2d waveguide with bounded inhomogeneity, b) poles of the Laplace transform of 
the different waveguide modes 



The first row in these matrices correspond to the boundary degree of freedom 
uq. The local clement matrix for the infinite clement is then given by 

Note that it is linear in k 2 and the matrices are explicitly known. In the 
equivalence of the variational equation in iJ 1 ([0,a]) x H + (D) and the classical 
problem ^ is shown. Moreover, the stability of the Hardy space infinite element 
method follows with a Garding inequality and exponential convergence in the 
number of degrees of freedom TV for the Hardy space is proven. 

Remark 4. In the space domain the monomial basis functions z 3 correspond 
to the functions 

Therefore for the one dimensional problem an optimal choice for scattering 
problems is kq = k since in this case the exact transparent boundary condition 
is obtained even with no degrees of freedom in H + (D). For resonance problems, 
kq should be chosen in the region of the complex plane where resonances are of 
interest. 



2.2. Tensor product elements for a semi-infinite strip 

For the multi-dimensional case we first present the extension of the one- 
dimensional Hardy space method to Helmholtz problems for waveguides with 
locally bounded inhomogeneities like the circle in Fig. 2(a) Given some incom- 
ing wave ui satisfying the homogeneous Helmholtz equation — Au; — n 2 Ui = 0, 
the total wave u is the sum of u; and a scattered wave u s , which has to satisfy 
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a radiating condition. The problem is given by 

- Au(x,y) - n 2 u(x, y) = 0, (x,y) € O C R X [0,7r], 
<9„it = on 90, 
M K0 £u 8 (.,y) e #+(£), 2/e[0,7r]. 



(12a) 
(12b) 
(12c) 



The Laplace transform in (12c I is applied for fixed tangential variable y to the 
scattered wave u s outside the gray shaded domain, i.e. for x in both directions 
±oo. Note, that u s can be written as a superposition of one-dimensional waves 
with wavenumbers k n := V k 2 — n 2 



u s (x,y) 



n=0 



i cos(ny)e lknX . 



(13) 



As sketched in Fig. 



2(b) there exist only a few guided modes and infinitely 
The parameter kq in the Mobius transform has to 



many evanescent modes, 
be adapted to the region, where the Laplace transforms for all these modes 
are holomorphic. This way the results of [T3l [11] carry over to the waveguide 
problem (12), even though this is not directly contained in these papers. 

fa, 6] x [0, tt] of Fig. 



In the gray shaded domain Q- mt :— f2 (~l 



we use a 

standard finite element method to approximate the total wave u. Since only the 
scattered wave satisfies the radiation condition, we have to solve the problem 
in the exterior domains for u s and not for u. On the artificial boundary T, that 
separates Vt- mt from U, cxt this gives rise to a jump condition for the Dirichlet 
values as well as an additional boundary term from the jump in the Neumann 
values: 



(Vzt • Vu — n 2 uv) d(x, y) + / (Vu s • Vu — k 2 u s u) dxdy 

Oi n t J J — OO 

7T />00 



(Vm s ■ Vv - k u s v) dxdy = / {{d x u{) (b, y)v(b, y) - (d x Ui) (a, y)v(a, y)) dy 

Jb JO 

for suitable test functions v. The infinite integrals can be transformed into the 
Hardy space using ^ and as in the one-dimensional case the decomposition 
([6|, which ensures the continuity of the solution over the interfaces. E. g. for 
x > b the stiffness and mass integral become 



Vu s ■ Vv dxdy = -2m a I A[T+[ ) ,T+ ( J?^\ ) j dy 



7T roc 



U(;y)J" + \V(;y) 



K Jb 



\U(;y))> '- \V(;y) 



Since on the interface [0, it] there already exists a discretization consisting of the 
traces bm^ of finite element basis functions in the interior domain, we use these 



8 



basis functions for the boundary values u s0 and vq- 



u s0 (y) = c m 6 m(y). v e [o,tt]. 

m=0 



For [/, V € _ff + (D)(8iiJ 1 / 2 ([0, 7r]) it is reasonable to use tensor product elements: 
U(z,y)=J2Y / c m , i z i b^(y), z&S\ y€[0,n]. 

m=0 i=0 



In this way the integrals ( 14 ) give rise to tensor products of the boundary 
matrices 

= f 9,6^) (y) c^ 5 (y) d», M bd = f ^ (y) b& (y) dy (15) 
Jo Jo 

and the Hardy space matrices 5 HSM = -2iK T^ >+ T N ,+ and Af HSM = ^T^_T N - 

S cxt = S mM ® M bd + M HSM <g> S bd , M ext = M HSM <g> M bd . (16) 

£.5. Tensor product elements 

In the general multi-dimensional case we consider the Hclmholtz equation 

- &u(x,y) - K 2 n(x,y)u(x,y) = 0, (x,y)eCl, (17) 

with a potential n, an unbounded domain SI, with some boundary condition on 
dfl and as a radiation condition the pole condition along a generalized radial 
direction for u. Soon it will become clear what is meant by a generalized radial 
direction of u The assumptions on the potential n will be given in Remark [7] 

In 1 11 j the computational domain is obtained by intersecting Jl with a ball 
B a of radius a, f2i nt := B a nf2, such that the unbounded exterior is f2 cx t := R d \ 
S a . Using polar coordinates in f2 ext and separation of variables the unbounded 
radial direction and the bounded surface directions separate. Hence, the one- 
dimensional approach can be applied to the radial part of the exterior solution 
and a standard finite element method handles the interior part as well as the 
bounded surface directions. 

However the boundary of f2 cx t need not be a sphere, arbitrary convex poly- 
gons P can be used to split Q, into fi e xt := O \ P and f2; n t :— P n VI with 
interface Y := dP. Here for simplicity we present only the two-dimensional 
case. In Fig. |2.3| two different segmentations of the exterior domain are illus- 
trated: The left one decomposes the exterior into semi-infinite strips and infinite 
triangles, while the right one uses semi-infinite trapezoids. 
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Figure 3: different segmentations of the exterior domain 



2.3.1. Infinite triangles and strips 

We first present the Hardy space infinite element method for the segmenta- 
tion using infinite triangles and strips. For the semi-infinite strips we use the 



elements from section 2.2 Hence only the implementation of the method for 
the infinite triangles is considered below. If P is the vertex of such an infinite 
triangle and n\ and are the unit normal vectors of the neighboring strips, 
then the triangle is given by T = g ([0, oo) x [0, oo)) with the linear mapping 

S^,v) = P + ^n 1 +r)n 2 , G [0,oo) x [0,oo). (18) 

If we define ^(£,77) := u(g{^,rj)) 011 the reference, then with the constant Ja- 
cobi matrix J of g, the mass and stiffness integrals for the infinite triangle are 
transformed according to 



u vd(x,y) = iiv \J\d(£,ri), 

T J[0 : oo)x[0,oo) , lg . 

/ V xy u ■ V xy vd(x,y) = / J' T \7^ v u ■ J~ T \7 (v v\J\ d(£,r)). 

JT J [0,oo) x [0,oo) 

In contrast to the strips, the integrals include the Jacobi matrix and there are 
two infinite directions to which ^ is applied to. If we define the constant matrix 
G := |J|J _1 J _T , then the local element matrices for each infinite triangle T 
are given by 

S T =G n S HSM ® S mM + G 12 S mM ® M HSM 

+ G 21 M HSM (8 S HSM + G 22 M HSM ® M HSM , (20) 
M T =|J|Af HSM ®A/ HSM . 

Remark 5. In the transformation g we use the unit normal vectors, which 
guarantees the continuity of the solution along the infinite rays. 
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(0.1 1' 
4 



(0,0) (1,0) rr 

1 2 ' 



Figure 4: Transformation of each trapezoid 



2.3.2. Infinite trapezoids 

The method presented in |2.3.1| has the advantage to be easy to implement, 
but it results in extra degrees of freedom in the infinite triangles. The second 
methods avoids these degrees of freedom by using infinite trapezoids. The infi- 
nite rays, which are no longer normal to the boundary, could e.g. be constructed 
in 2d with bisecting lines. Another possibility is to choose a reference point Pq 
in the interior domain and construct the rays R for each vertex V of the bound- 
ary by R = V — Pq. General conditions for suitable segmentations in infinite 
trapezoids may be found in [531 I2H US] • 

Given a segmentation with finite trapezoids, the trapezoid element is the 
image of a reference rectangle (see Fig. [4]) under the affine bilinear mapping g 
with 

(x,y) = g(?7,£) = R o Q(r],0 + (zi, Vif , 



(21) 



where 



(x,y) = Q(t],0 



b£, + (a 



(22) 



with h v = v /(x 2 - xi) 2 + (y 2 - yi) 2 , 
y3) T /||(a:4 ^ x 3 ,y 4 - g/3 ) 1 1 2 , b = (x 3 - x 4 ,V3 
X3,V4 - 2/3) H2 and h{ = y/\\(x 3 
signed distance variables. 



a = {x4 - x 3 ,y 4 - y 3 ){x2 - x 3 ,y 2 
Vi){x i - x 4 ,yi - y 4 ) T 



X2,V3 - V2, 



X4 — 

Note, that a and b are 



Remark 6. The variable £ plays the role of a generalized radial variable, whereas 
77 is the surface variable on T. To guarantee continuity of the discrete solution 
in the exterior domain it is important, that the radial variable £ along the rays 
of the segmentation is independent of the neighboring infinite elements. This is 
fulfilled, if two neighboring (finite) trapezoids have the same boundary vertices. 



The rotation R is given by 
(x,y) = R(x,y) = 



X 2 
2/2 



x\ yx 
y\ x 2 



in 



^(x 2 - xi) 2 + (y 2 - yi) 2 
The Jacobi matrix J of the transformation g and its determinant are 



(23) 



J 



hr, + (a + b)£ -b + (a + b)r) 
h£ 



\J\=h 6 (h v + (a + b)£) (24) 
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and the inverse of J is 



J- 1 



1 



h v + €(a + b) hz{h v + £(a + b)) 



b-(a + b)i] \ 
) 



(25) 



Note, that J is no longer constant. Mass ans stiffness integral transform as 

u v\J\ d(rj,£), 



in (191 



u v d(x, y) 
d(x,y) 



[0,l]x[0,oo] 



/ J- T V^u- J- T V^v \J\d(r),£). 

J[0,l]x[0,ool 



(26) 



As before we use tensor product basis functions on the reference element & m j := 
&m <8>&i . Again the integrals decouple, such that for the mass integral we obtain 

/ b m>i b nd d(x, y) = / dr} / bf^bfh^hr, + (a + 6)0 d£. (27) 

JT JO JO 

Due to the inverse Jacobian J -1 the stiffness integral is more complicated. We 
obtain 

/ ^ xybm,i ' ^ xybn,j d\X, y) — 
JT 

d v b$ (^ + (6-(a + 6)7 7 ) 2 )9,SWdry 



bfbf 



o ^(^ + (0 + 6)0 



d£ 



For functions depending on ?y we use again the traces of finite element function 
in f2; n t. Hence, the bounded integrals over rj can be treated in the usual way by 
quadrature formulas. For the infinite integrals we apply the identity ^ and the 
decomposition ^ to transform the ^-direction into the Hardy space H + (D). 

We have to take care of the factors £ and (£ + c) _1 with a constant c > in 
(27) and (28). For these we need one additional operator V : H + (D) — > H + (D) 
implicitly defined by 

M K „C{(')f} = M K0 {- (£/)'} = V (M K0 Cf) . 
Direct calculations yield 

(VF) (z) = { A F ^F'(z) + Z —±F(z), F e H+{D). (29) 
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space domain 


M K0 C 


def. 


implementation 


/ 


—T-(f ,F) T 


© 


(11 


1 


/' 


r + (/ ,F) T 






( 


11 


) 


(•)/ 


V{M K0 Cf) 


(29 


) 


( 


30 


) 


• + c 


(V + c idy 1 (M Ka Cf) 


(31 


) 


num. inv. of {T>n + c id) 



Table 1: Operators for the HSIE method; for the basic identity see |7]) 



If we use the set of trigonometric monomials up to the order as basis functions 
in H + (D), we get the discrete operator 



V 



N f 



( -1 1 

1 -3 2 

2-5 3 



V 



-2N £ - 1 J 



(30) 



Obviously it holds for the operator (T> + c id) 1 
M,.,C ' ' 



/ } = (V + c id)- 1 (M K0 Cf) 



c> 0. 



(31) 



Note that both operators and the matrices "Dn 6 and (2?at £ + c id(;v£+i)x(iVe+i)) 
are symmetric. 

With the operators given in Table [I] and the identity |7]) we are able to 
transform all integrals over £ into the Hardy space H + {D). Using a Galerkin 
ansatz there with monomial basis functions, local matrices are obtained. For 



example the integral in (27) yields the local mass matrix 
2hf 



M := 



^T]f (h v id(7v (:+ i)x(Ar 5 +i) + (a + b)V Ni ) Tjv c ,-. 



IK 



(32) 



In the same way we can treat the integrals in ( 28 ) and obtain the matrices 



Lqo 


•- . , T N _ (h v id(Ar ?+ i; 


x(iv e +i) 


+ (a + b)V N( ) 1 T N(1 - 


Lqi 


:= -2 T Nt _T N(t+ , 






Lio 


'■= -2 T N(t+ T Nf -, 






in 


— h 1 N e ,+ IN id(JVe+i)> 


(W 5 + l) + 


(a + b)T> Ni ) T N(t+ . 



(33) 



Remark 7. The method is applicable to a wide range of scalar Hclmholtz- 
type problems. Since each segment in fi e xt is treated separately, unbounded 
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inhomogeneities as e.g. waveguides are possible. Even coefficient functions n 
with unbounded support are possible, if there exist a segmentation of 51 cxt such 
that in each segment the function can be written in sums of terms (£ + a j) kj 
(including negative powers) and functions depending on 77 

n(£,v)=no g (£,r)) = ^2{£ + a j pc j {r)), a s > for k 3 < 0. (34) 

Remark 8. The exact statement of the tensor product space is due to the de- 
composition in ^ a little bit complicated. Another problem is, that the infinite 
integrals have to be bounded. This can be done by choosing test functions v, 
which decay fast enough and who are dense in the Hardy space after transfor- 
mation. The details can be found in [XT] . 



3. Perfectly Matched Layer 

Exterior complex scaling was introduced by Simon [53] to facilitate the math- 
ematical formulation of boundary conditions for the wave functions in quantum 
mechanics. It is shown by Chew and Weedon [3] that Berenger's [2j PML devel- 
oped for transient Maxwell's equations may be interpreted as a complex scaling 
of the exterior solution. Thus PML can be regarded as equivalent to exterior 
complex scaling. 

Starting from a discretization of the exterior as in Fig. |8(b)| the generalized 
radial coordinate is scaled by a constant complex factor a, such that scattered 
outward radiating waves are damped exponentially. Using this discretization 
no special corner conditions are required. The exponential damping justifies 
to truncate the infinite domain some distance away from the boundary of the 
computational domain f2j„ t and to impose homogeneous Dirichlet or Neumann 
boundary condition. This way the computational domain is surrounded by a 
finite layer. The convergence of the PML method for homogeneous exterior 
domains is analyzed by Lassas and Sommersalo [18] for the scattering problem 
and by Kim and Pasciak [17] for the resonance problem. 

In our numerical experiments we compare the HSIE method with the adap- 
tive PML method described in [27] [2J]- The special feature of this PML is that 
the thickness of the layer is chosen adaptively based on an a posteriori estimate 
of the error introduced by truncating the layer and taking into account the dis- 
cretization error of the interior. The distribution of the grid points is based on 
the observation that inside the PML short waves that require a fine grid to be 
resolved with a certain accuracy are damped much faster than long waves. Long 
waves in turn are well resolved on rather coarse grids. 



4. Convergence test: A strip waveguide 

In this section we compare different discretizations of the exterior domain 



for the Hardy space infinite element method applied to (17) with k 



2tt 
1.5' 



the 
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Figure 5: a) real part of n;; b,c) 1st mesh with two different exterior discretizations; d) 2nd 
mesh 



refraction index 



n(x,y) = n(y) 



n\ y£ (-a, a) 
n\ y€R\(-a,a) 



(35) 



and a = 0.0365, ri\ = 1.45 and n 2 = 3.4. The incoming wave is given by 
Ui(x,y) — v(y)e lli * x , with 



v(y) 



' de^-"'" 2 " ,y<-a 
C 2 e- l V n l K2 - K ly + Cse^^-^y , y € (-a, a) (36) 



for k x > and complex coefficients C\, C2, C3 and C4, which have to ensure 
the continuity of v and v' in HL 

Remark 9. The incoming wave should satisfy the Helmholtz equation. Plug- 



ging the ansatz m(x,y) = v(y)e lK * x into (17 1 leads to the eigenvalue problem 



-d,, — K^n) v 



-kIv 



for the eigenpair (k%,,v) elx H 2 (M.). If the jump in the refraction index n is 
large enough, such an n x £ (riiK, n 2 K) exist and the corresponding eigenfunction 
v is exponentially decaying for y — ¥ ±00 and oscillating in (—a, a). 

We solved the problem for strip waveguide for the two different meshes in 
Fig. [5] The incoming wave is coupled via the jump conditions described in 
section |2.2| on the left vertical boundary part of the domain, which is chosen 
sufficiently large, so that u\ can be set to in the lower and upper left corner. 
In order to test the Hardy space method the right vertical boundary is for one 
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Figure 6: Relative H 1 ^^) error of the Hardy space method vs. the number of degrees of 
freedom in radial direction 



part of the computations very small. In the area around the right waveguide 
port the boundary (and with it the Hardy space method) has a big influence on 
the numerical solution, which should approximate the incoming field U{. 

For the exterior domain two different types of discretization are used: First 
we combined infinite strips with infinite triangles (Sec. |2.3.1 1. Second we used 
the trapezoids of section 2.3.2 As shown in Fig. [5] the rays are chosen in 



normal direction to the boundary and in corners as bisecting lines. For the 
computations we refined the diagrammed meshes once and used a high order 
method with polynomial order 7 for the finite element method in the interior 
domain. 

Fig. [6] shows for the two different discretizations and the two meshes in 
Fig. [5] exponential convergence of the Hardy space method with respect to the 
number of degrees of freedom in radial direction. For the rectangular mesh 
both discretizations give the same result, since the influence of the upper and 
lower right corners on the discrete solution is very small. For the challenging 
mesh the trapezoidal discretization needs approximately one degree of freedom 
more than the non-uniform discretization. The computational costs for both 
discretizations are similar. 

Fig. [7] shows the logarithm of the absolute value of the degrees of freedom 
in the Hardy space. The white domain in the middle is the interior domain and 
we have plotted only the degrees of freedom for the trapezoidal discretization. 
As expected from the one-dimensional case and from the theory about spherical 
exterior domains they decrease exponentially. Nevertheless, they do not decrease 
uniformly in all directions. Especially near to the exit port of the waveguide the 
decreasing factor is lower than in the other regions. For this reason a strategy to 
choose the number of Hardy modes adaptively is currently under investigation. 
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Figure 7: Logarithm of the absolute value of the degrees of freedom along the rays of the 
exterior domain for two different meshes in Fig. [5] 



5. Numerical example: a micro cavity resonator 

This example is taken from Hammer [S] and consists of two waveguides 
coupled by a cavity. In order to exclude discretization effects originating from 
the resolution of the layout a geometry is chosen that is defined by polygons. 
The computational domain has a size of 3.5/im x 4.546/im, with a square cavity 
a = b = 1.451/mi in the center. We have c = 0.2745^m and d — 0.073/xm in 



Fig. 8(a) 



5.1. Scattering problem 

To model the scattering of an incoming wave by an object the incoming 
waveguide mode of the last section is coupled by a jump condition of the Neu- 
mann and Dirichlet data at the left vertical boundary part for two different 
wavenumbers k\ = jfe and fc 2 = , .-i 7 ^ — . fc 2 is close to a resonance, whereas 

1 1.5/jm ^ 1.5759um ^ ' 

for ki the cavity has only little effect, which can be seen in Fig. [9J For fc 2 the 
wave propagates through the cavity into the lower waveguide. 

For these computations we used for the interior finite elements of 5th order 
and refined the coarse grid |8(b)] three times uniformly. In the exterior domain we 
used the trapezoidal Hardy space infinite element method with the parameter 
«o = 8 + 5i and 30 degrees of freedom in order to discretize the Hardy space. 
In total we got approximately 200.000 degrees of freedom. 

5.2. Resonance problem 

Starting from a coarse mesh the interior is adaptively refined using a resid- 
ual based error estimator [27]. On each refinement level the eigenvalue of the 
resonance problem is calculated. To evaluate the relative error in the eigen- 
value, a reference resonance frequency if uj = 1.1951173e + 15 — 1.489202e + 13i 
is calculated on a very fine mesh using the PML. The wavelength of the reso- 
nance is A = 2ttc/uj, where c = 299792458m/s is t he speed of light. Hence the 



wavenumber ki = 2n/X-i = . — of Section 5.1 is close to the resonance. 
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Figure 8: a) Schematic geometry of the micro cavity resonator, b) Coarse grid discretization 
by triangles in the interior and trapezoids in the exterior, (a = b = 1.451/im, c = 0.2745/im 
and d = 0.073/im) 
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Figure 10: Comparison: PML (dashed line with filled markers) with HSIE (solid line) with 
Krj = 5 + 3i for various finite element degrees: 1 (o), 2 (A), 3 (□), and 4 (0). Left: Work- 
precision diagram showing the relative error in the eigenvalue vs. cpu-time in seconds. Right: 
Convergence of PML and HSIE showing the relative error in the eigenvalue vs. number of 
degrees of freedom. 



GSMBSDOOO OOOOOOOOOOO 
10" 4 ^ 



10"° - 

iy<xw vv vvvvvvvvvvv 
io- 8 <4 



10 

10" 1 



\DD □□□□□□□□□□□ 

00 OOOOOOOOOOO 



10 20 30 40 50 

N coefficients 



Figure 11: Rel. error in the eigenvalue vs. number of Hardy-modes (ren = 5 + 3i) for finite 
elements of degree p = 1 (o), p = 2 (A), p = 3 (□), and p = 4 (<))■ 



In this example two methods, PML and HSIE, to realize transparent bound- 
ary conditions are compared. The results for the PML are obtained using the 
adaptive PML as described in Section [3] The results for the HSIE are obtained 
using ./V = 2, . . . , 50 Hardy modes and selecting the best result. [10] shows the 
error versus the cpu time (left) and the total number of degrees of freedom 
(right) that is required to solve the eigenvalue problem on the finest refinement 
level. The initial guess is 1.195 • 10 15 - 0.01489 • 10 15 i. The HSIE method in 
this example yields results that are better or at least as good the results ob- 
tained with the PML method. The implementation was done in the C++ code 
JCMsuite [IS]. 

Fig |ll| shows the super algebraic convergence in the number of Hardy modes. 
For these calculations kq = 5+3i. The number of degrees of freedom to discretize 
the interior is about 944000 for p = 1, 169000 for p = 2, 121000 for p = 3 and 
103000 for p = 4. For low accuracies very few Hardy modes are required, 
e.g. to reduce the approximation error of the transparent boundary condition 
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below 10 -6 7 Hardy modes are sufficient. To obtain a relative error of about 
fCP 10 almost one third of the total number of degrees of freedom is spent on 
approximating the transparent boundary condition. 

6. Conclusions 

We have presented here the Hardy space infinite element method from a 
practical point of view. In this form inhomogeneous exterior domains can be 
treated as well as arbitrary convex polygons as artificial boundaries. The nu- 
merical results show superalgebraic convergence with respect to the degrees of 
freedom in the Hardy space, i.e. the degrees of freedom in the radial direc- 
tion. The method preserves the eigenvalue structure and is therefore well suited 
for solving resonance problems. Compared to former realizations of the pole 
condition like the cut function approach f }25|) it is not possible to recover the 
solution in the exterior domain directly. In this point the HSIE method behaves 
like perfectly matched layer methods: The degrees of freedom in the exterior 
domain are physically irrelevant. They just provide for a good approximation 
at the exact transparent boundary condition on the artificial boundary. 

A comparison of the results of the HSIE and the used PML is already given 
in the last section. From the point of implementation the HSIE requires a new 
(in) finite element, whereas the PML only changes the variational formulation 
of the problem. On the other hand the HSIEs show exponentially convergence, 
while the PML inherits the convergence order of the used finite element method. 
One basic difference between both methods is given by the nature of the dis- 
cretizations. For the PML there exist two steps: First truncating the infinite 
PML-domain and than discretizing the finite layer using finite elements. The 
HSIE uses a (transformed) variational formulation of the whole infinite domain. 
The only discretization results from the Galerkin method. 

The HSIE is not restricted to scalar, time-harmonic problems. In [5T] a sim- 
ilar version of the method is used for solving time-depending problems. More- 
over, there exist first results for Maxwell's equations. 
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